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Abstract 

We report findings related to a two dimensional viscous fingering problem solved with a timespace method and 
anisotropic elements. Timespace methods have attracted interest for solution of time dependent partial differential 
equations due to the implications of parallelism in the temporal dimension, but there are also attractive features in the 
context of anisotropic mesh adaptation; not only are heuristics and interpolation errors avoided, but slanted elements 
in timespace also correspond to long and accurate timesteps, i.e. the anisotropy in timespace can be exploited. We 
show that our timespace method is restricted by a minimum timestep size, which is due to the growth of numerical 
perturbations. The lower bound on the timestep is, however, quite high, which is indicative that the number of 
timesteps can be reduced with several orders of magnitude for practical applications. 
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1. Introduction 

Many engineering problems are governed by partial differential equations which can only be solved numerically. 
The prospect of parallelising such simulations in the time dimension have attracted attention by the scientific commu¬ 
nity dEHIl, and the parareal algorithm la is popular and widely studied. The parallelism in the time dimension 
can be achieved by leaving the distinction between time and space, in what is normally referred to as timespace 
methods. The technique is often used for solving hyperbolic problems with dicontinuous Galerkin (DG) formulations 
i), as demonstrated for the full four dimensional case in the context of rotor flows m For that application all the 
complications of a time varying geometry disappear, when switching to a timespace method. This advantage was also 
exploited for dynamic stall prediction as simulated with the Navier-Stokes equations M- That work utilised spacetime 
adaptation and, as it is often the case, a compromise was made in the form of timeslabs, which allows for a tunable 
problem size (91. 

Anisotropic mesh adaptation is an established technique for ensuring computational efficiency in the context of 
multiscale problems oniiiiiiniini. There exists many heuristics methods for applying the technique to transient 
problems iniiiiEi, but the only way to exploit the anisotropy in the physics as seen in timespace, is by using a 
timespace method. This is also a conceptually simple way to avoid the interpolation errors and heuristics associated 
with combining conventional timestepping algorithms with mesh adaptation. The benefits of anisotropic mesh adap¬ 
tation generally increases with dimensionality due to the fact that elements can be stretched in more dimensions, and 
one could thus expect significant benefits from the use of a four dimensional anisotropic mesh adaptation algorithm, 
but (to our knowledge) anisotropic mesh adaptation has only been demonstrated for two and three dimensions. We are 
not aware of any studies involving the use of three dimensional anisotropic meshes for solving a timespace problem. 

In this work, we focus on a miscible viscous fingering problem, which is transient and has two spatial dimensions. 
The problem is known to exhibit chaotic behaviour ca, which causes some issues for the combination of a timespace 
method and anisotropic mesh adaptation. We attribute this to the fact that the numerical noise can grow in time. This 
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could mean that a non-linear solver is unable to converge to the error level that one would expect for the computational 
resources utilised. We suggest to use the method of timeslabs, where the slab thickness is a kind of timestep size. The 
method is used to investigate the relation between the timeslab thickness, the computational resources used and growth 
of numerical errors. 

2. Anisotropic Mesh Adaptation 

Mesh adaptation is the art of choosing the discretisation that minimises the error level for a given computational 
cost. It is important to distinguish between p- and h-adaptation, which relate to the discretisation order and length 
scale, respectively. Anisotropic mesh adaptation is the most general type of h-adaptation, because it not only considers 
the local element size, but also the element orientation. The optimal size and orientation can be expressed in terms of 
a metric tensor field [El, At . This has units of inverse square length, and it can be used to map elements to metric 
space, where the optimal element has unit edge lengths. The metric that minimises the interpolation error iflTl ofa 
variable with Hessian, H, is 

^ ^ (det[ate(H)] j ate(H), (1) 

where d is the number of dimensions, q is the error norm to be minimised, det is the determinant, cr is a scaling factor 
and ^ takes the absolute value in the principal frame. The metrics of several variables can be combined with the 
inner ellipsoid method ifTSl . illustrated for two dimensions in figure 



Figure 1: The inner ellipse method is illustrated in the case of intersection, but it is common to see one ellipse entirely within the other, such that 
anisotropy is preserved. 

We use 1st order polynomials for all fields, so we compute the Hessian by first performing a Galerkin projection 
of the gradient onto 1st order polynomials and then repeating the process to get a nodal Hessian. The nodal metric 
can then be calculated explicitly using equation Q- 

We use the technique of local mesh modifications to adapt the mesh to the metric. There are four operation types: 
Coarsening, swapping, refinement and smoothing, as shown in figure Coarsening removes edges that are short in 
metric space without regard to mesh quality, but the other modifications only allow the worst local element quality to 
go up. The mesh quality is quantified by means of the Vassilevski functional iTT^ . 


Coarsening Edge swapping Refinement Smoothing 



Figure 2: Four local mesh modifications are illustrated: Coarsening, refinement, swapping and smoothing. Only the edge swapping operation is 
illustrated in three dimensions, as it is trivial to generalise the other operations to three dimensions. Coarsening is the only operation that is allowed 
to reduce the worst local element quality and therefore edges marked for refinement are not guaranteed to be split. We do not use face to edge 
operations also known as 2-to-3 operations. 

We use an Octave/MATLAB implementation, which is suboptimal in terms of performance, especially compared 
to a similar C-f-f implementation CD. 
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3. Setup 


We consider a simple two-phase miscible problem (161, where the viscosity, rj, depends exponentially on the 
saturation, 0, with (p = 0 and 0=1 corresponding to the low and high viscosity fluid, respectively. 


T]((p) = 


( 2 ) 


where ^ determines the viscosity ratio. The velocity, v, is given by the Darcy equation. 


V = 


— V/7, 
riW ^ 


(3) 


where p is the pressure and k is the permeability. Mass conservation leads to a Poisson equation for the pressure. 


o = v-(pv) = 




(4) 


where p is the density. This drops out of the equations, because we choose to study the case of equal densities for the 
two fluids. The saturation is convected with the velocity given in equation <0, 


dd) 

0 = ^ + V . V0 

ot 

dd) K ^ ^ 

= ^ 
dt 


(5) 


The equation is treated as a convective equation with a three dimensional velocity vector, and it is stabilised with 
streamline upwind Petrov/Galerkin diffusion. The local characteristic length scale for the mesh is calculated by 
computing the Steiner ellipsoid for every element and projecting this along the velocity direction. 

The governing equations 0 and 0 are solved with the boundary conditions illustrated in flgurej^ The use of a 13- 
sided polygon to approximate a circle avoids extrapolation, when comparing solution fields on different meshes. The 
perfect circular initial condition for the saturation makes the solution undefined, which causes the growth of numerical 
perturbations to become very obvious. All equations are solved with the finite element method using FEniCS (20l, an 
open source finite element package. Both pressure and saturation is discretised with first order polynomials. A direct 
solver (LU) is used for equations Q and while an iterative solver (CG-fILU) is used for Galerkin projections. All 
computations are single threaded. 

We non-dimensionalise the equations with A/? as characteristic pressure and R as characteristic length scale. 


R^tjo 

X = R% t = -f, p = App, 

VcharKAp 

r = OAR, ri = 0.2R, r 2 = 0.3/?, ^ -2, ffinai = 1 


where x, i, p are dimensionless space, time and pressure, respectively, ffinai is the simulation time. Vchar is a numerical 
parameter, which determines the relative resolution of space and time. We fix it at unity. Another numerical parameters 
is the norm of the interpolation error to be minimised, and we go with the popular choice of ^ = 2 for this. 

The saturation at the initial time is imposed by means of a Dirichlet boundary condition. We initialize the non¬ 
linear solver with 0 = 0 and p = pinit, where 


r-0.1 
“ 0.9 ’ 

f being the radial coordinate. We use a segregated fixed point method to deal with the non-linearity of the equa¬ 
tions. The timespace problem is solved in slabs with a thickness of At (omitting the tilde), i.e. 

#0 Initialise mesh, and set 0 = 0, p = pinit- 
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Figure 3: The initial condition for the viscous fingering is sketched with the dark grey fluid pushing out on the light gray fluid, which has a higher 
viscosity. The initial saturation is unperturbed and it is thus primarily the length scale of the initial mesh that triggers the instability. Zero pressure 
is prescribed at the outer boundary, while p = Ap and 0 = 1 are prescribed at the inlet. 


#1 Solve equation Q for the pressure with fixed saturation. 

#2 Solve equation Q for the saturation by taking the viscosity to be constant using the previous saturation. 

#3 If it is the first iteration at this timeslab, go to #1. Otherwise, we calculate metric fields associated with the 
pressure and the saturation. These are combined with the inner ellipsoid method (see figure [^, the mesh is 
updated and the fields are interpolated onto the new mesh. If it is not the first timeslab, the mesh is fixed at the 
boundary matching up to the previous timeslab, 

#4 If 20 iterations have been carried out with the current timeslab, we proceed to the next one by mirroring the 
mesh in a plane normal to the time direction, and set 0 = 0, p = pinit- We then go to #1. 

The mesh is changed in every iteration, and therefore the non-linear solver cannot converge to machine precision. A 
typical timeslab is shown in figure [^together with the result of a full timespace simulation. Both figures illustrate the 
value of full timespace anisotropy. 


4. Results 

We show result for timeslab thicknesses At = 0.05, 0.1, 0.2 and 0.5 with cr = 0.01 as well as for mesh tolerances 
cr = 0.04, 0.02 and 0.01 with At = 0.1. The result at the end time is plotted in figures and respectively. The 
final iteration number, time. At value, cr value and final total computational time is printed above each image. As one 
would expect, the shape and number of fingers clearly varies with cr, but not with At. 

We measure convergence in terms of the difference between saturation fields in consecutive iterations, 

2 ^ 

' " ■ 

This is plotted in figure [7Ja-b). Figure |7Ja) clearly shows that simulations become more accurate as At is lowered, 
but the convergence is different when cr is decreased, at least in the last half of the simulations: The error seems to 
scale with cr, but it scales down to a certain level with At. By comparing figure [TJa) and (b), we can see that this is the 
level of the overall discretisation error. We attribute this effect to the fact that perturbations due to the mesh grow with 
time due to the non-linear nature of the problem. The extreme case of this is illustrated in figure where snapshots 
of the 0 = f = 0.5 contour is shown for three different iterations with At = 0.5. There are significant differences due 
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Figure 4: A simulation with cr = 0.02 and At = 0.1 is plotted in terms of a wireframe of the final iteration with slab #5 (left) and all slabs in terms 
of the 0 = 0.5 isosurface (right). The time dimension is from left to right pointing out of the plane. Large elements clearly show up away from the 
interface, but even at the interface the anisotropy in timespace can be exploited such that very few elements are needed. Note that the mesh is fully 
tetrahedral; the quadrilaterals on the isosurface representation appear, when cutting tetrahedrons with two nodes on both sides of the threshold. 


to the growth of mesh perturbations, particularly between iteration 16 and 20. This is to say, that the solution is not 
well defined and therefore it has an oscillatory component arising from the changing mesh. The error associated with 
this component increases with time, so it can be smaller than the overall discretisation error, if At is sufficiently small. 
Focusing on this case with cr = 0.01 and ^ = 2, it looks like Af = 0.1 is a good value. 

5. Conclusion 

The results show that it is possible to simulate chaotic phenomena with a timespace method. A CFL like condition 
still seems to exists, so the timeslab cannot be too thick. This has to do with the fact that the chaotic nature of the 
problem allows for a range of solutions and the difference between these solutions, should not exceed the overall 
discretisation error. The lower limit for the timeslab thickness, is high enough to indicate that the total number of 
timesteps can be reduced by orders of magnitude compared to a conventional method. 

6. Suggestions for Future Work 

It would be interesting to compare an established method against an optimised timespace implementation in the 
context of large scale simulations and anisotropic mesh adaptation. If practical problems are to be solved, a four 
dimensional anisotropic mesh adaptation algorithm will ultimately be needed. 
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Figure 5: The saturation at the final time is plotted for four simulations with cr = 0.01 and At= 0.05, 0.1, 0.2 and 0.05. The shape of the fingers 
vary, but the number and length seems independent of At. 



i = 200, t = 1.0, a = 0.01, At = 0.10, 14.3h 
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Figure 7: The square root of the residual is plotted as a function of normalised iterations numbers, (a) shows varying Af at cr = 0.01, while (b) 
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